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Crystal Nucleation of Colloidal Suspensions under Shear 
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We use Brownian Dynamics simulations in combination with the umbrella sampling technique to 
study the effect of shear flow on homogeneous crystal nucleation. We find that a homogeneous shear 
rate leads to a significant suppression of the crystal nucleation rate and to an increase of the size of 
the critical nucleus. A simple, phenomenological extension of classical nucleation theory accounts 
for these observations. The orientation of the crystal nucleus is tilted with respect to the shear 
direction. 
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The formation of crystals in a supercooled melt is a fas- 
cinating yet complex process. It is initiated by a micro- 
scopic nucleation event. The resulting embryonic crys- 
tal then grows to macroscopic size. Understanding the 
principles of nucleation and growth is essential for many 
applications rangingfrom tailored protein crystallization 
to metallurgy 0,Sl3- ^^ present, the most detailed ex- 
perimental information on crystal nucleation comes from 
hard sphere colloids |j, la, IS, tD|- Such suspensions are 
ideal to study crystal formation, as the equilibrium and 
transportproperties of hard-sphere colloids are well un- 
derstood [Sj. Moreover, recent progress in computer sim- 
ulations has made it possible to predict the absolute rate 
of crystal nucleation in colloidal suspensions |3, ll3 and 
thus to compare with experiment. 

In the present Letter we explore the influence of shear 
flow on colloidal crystal nucleation. Note that applying 
shear is qualitatively different from the effect of pressure, 
temperature or additives, as the latter affect the thermo- 
dynamic driving force for crystallization or the rate of 
crystal growth. In contrast, a system under shear ends 
up in a non-equilibrium steady state. Several experimen- 
tal studies of the effect of shear on crystallization have 
been reported in the literature. Some of these report 
a shear-induced orderiiigof the liquid which enhances 
the nucleation rate HI [H El El , while others [HEl 
report the observation of shear-induced suppression of 
crystallization. Both phenomena can be qualitatively 
understood: on the one hand, shear may induce layer- 
ing in the meta-stable fluid, thus facilitating crystal nu- 
cleation. On the other hand, shear can remove matter 
from small crystallites and thus works against the birth 
of crystals. At present, it is not clear which mechanism 
is dominant, and under what conditions. In this Letter 
we combine the umbrella sampling technique from equi- 
librium Monte Carlo simulations with Brownian Dynam- 
ics simulations to study this non-equilibrium problem. 
We confirm that shear suppresses crystal nucleation, at 



least for small shear rates, as found by Butler and Har- 
rowell 17|, and in addition characterize the associated 
critical nucleus. 

Below, we consider homogeneous crystal nucleation in 
a simple model for charge-stabilized colloidal suspensions 
subjected to linear shear flow. The charged colloidal par- 
ticles interact via a repulsive Yukawa potential |S| 
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where n is the inverse screening length and r the mutual 
distance. The dimensionless strength of the interaction 
/3e has been fixed at a value /3e — 1.48 x 10"*, where 
(3 = \/{kBT) the inverse thermal energy and we used a 
cut-off at a distance 10/k. To model the time evolution 
of the sheared suspension, we used Brownian Dynam- 
ics jlSi Il9| . In this approach, hydrodynamic interactions 
between the colloids are ignored. This is justified at low 
volume fractions of charged suspensions. 

The Brownian-Dynamics equations of motion for a sys- 
tem in the presence of a steady shear rate 7 are of the 
form 

r,(i + 5t) = n{t) + St^ +Sr^ + Stjy,{t)x . (2) 

Here ri{t) = (xi{t) , yi{t) , Zi{t)) is the position of the ith 
colloidal particle at time t. In a small time interval St this 
particle moves under influence of the sum of the conser- 
vative forces fi {t) arising from the pair interaction (^ of 
particle i with the neighboring particles. During this mo- 
tion, the solvent exerts a friction. The friction constant ^ 
with the solvent is related to the diffusion constant D by 
£, — ksT/D, while the stochastic displacements are inde- 
pendently draw from a Gaussian distribution with zero 
mean and variance {{Sr^^)'^) = 2DSt, where a stands 
for one of the Cartesian components. The last term in 
Eq. Q represents the applied shear in the x-direction. 



and imposes an explicit linear flow field. For the simu- 
lations we used a cubic simulation box with 3375 parti- 
cles and Lees-Edwards periodic boundary conditions [20| ■ 
The total simulation time was up to 10'*/(k^Z?) for gath- 
ering statistics. The osmotic pressure P is kept at a 
constant value with isotropic volume moves. In practice 
this means that after a number of Brownian dynamics 
time steps the volume of the simulation box is attempted 
to be modified and the particles locations are scaled ac- 
cordingly. The resulting difference in potential energy is 
used either to accept the new volume or to reject and 
restore the old volume and particle locations, following 
the rules as used in normal Monte Carlo simulation of 
the isobaric ensemble jlq . The results for the zero shear 
case show full agreement with those we obtained by equi- 
librium Monte Carlo simulations. 

The number of particles inside the nucleus is deter- 
mined with the aid of bond-orientational order param- 
eters |21j, which characterize the neighborhood of each 
particle. By selecting particles with a solid-like environ- 
ment that are in each others neighborhood, all particles 
that belong to a cluster are identified. 

According to the bulk phase diagram, the stable equi- 
librium system would be a face-centered-cubic crystalline 
phase |lO|,l23 ^^^ '^^^ parameters. The system under con- 
sideration, however, is supercooled. Hence it remains liq- 
uid, even though the solid is more stable because, unless 
the nucleation rates are huge |23|, the simulation time 
required to observe spontaneous crystallization is very 
much longer than the duration of a run. Due to fluctua- 
tions the liquid will continuously form and dissolve small 
nuclei. Yet, the steady state probability P{n) that a crit- 
ical crystal nucleus of n particles will form spontaneously 
is extremely small. In order to speed up this process and 
obtain better statistics on the cluster size distribution, 
we used the umbrella sampling technique [2J| . The basic 
assumption underlying its usage is that the probability to 
find the system with a given cluster size is a unique func- 
tion of the thermodynamic state of the system and of the 
shear rate. To compute the probability to find the system 
in an unlikely state (such as a critical nucleus), we bias 
the Brownian-dynamics sampling in favor of the states of 
interest. The actual biasing procedure is identical to the 
one used in (meta-stable) equilibrium studies of crystal 
nucleation [9| , and merely works as a mathematical trick 
to measure the ratio of the function P{n), we want to 
obtain, over a known and fixed probability Pbias{n). All 
trajectories that are generated follow a normal path and 
are truncated by the bias when they deviate too much 
from the preferred cluster size. Rather than generating a 
new configuration, the last configuration is restored from 
which a new path is grown. Note that the trick of using 
umbrella sampling in a dynamical simulation is generally 
applicable in equilibrium and non-equilibrium situations, 
is not restricted to Brownian Dynamics, and enables one 
to obtain information on rare events. 
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FIG. 1: Negative logarithm of the probability P{n) of finding 
a cluster of n solid- like particles normalized by P(l) for pres- 
sure PP/ tC' — 0.24 and different applied shear rates, from 
bottom to top 7/(k^D) = 0, 0.8 x 10"^, 1.6 x 10"^, and 
3.2 X 10~^. The insets show typical snapshots of critical nu- 
clei for the largest shear rate and the zero shear case. 

After correcting for the biasing function, the cluster 
size distribution function is obtained. In Fig. ^ the loga- 
rithm of the probability function P{n) is shown for three 
non- vanishing shear rates. 

In the case that no shear is applied, one can relate the 
probability of finding a cluster of given size to the Gibbs 
free energy. It is therefore tempting to interpret the prob- 
ability functions as shown in Fig.^in terms of nucleation 
barriers |25(. Strictly speaking this is not allowed, since 
this idea stems from equilibrium considerations, while in 
the present case we treat a non-equilibrium system. How- 
ever, application of statistical mechanics outside equilib- 
rium can be useful (see e.g, |2y| for an effective temper- 
ature in a sheared system) and it is a challenge to check 
whether and to what extent equilibrium concepts are ap- 
plicable. In our case we consider the negative logarithm 
of the cluster size distribution function as an effective free 
energy. 

Under this assumption a simple extension of classical 
nucleation theory can be made, which incorporates the 
shear rate. In classical nucleation theory the Gibbs free 
energy AG of a spherical nucleus of radius R is given by 
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On the one hand there is a gain in energy proportional to 
the volume of the nucleus due to the difference in chem- 
ical potential A/i between the solid with density ps and 
the liquid phase. On the other hand we have a loss in 
energy, since an interface between the solid nucleus and 
surrounding liquid needs to be formed, described by 75^ 
the interfacial free energy. 

It is reasonable to expect that for moderate shear rates 
the chemical potential difference A/x and interfacial free 



energy ^sl will not be affected much. This would justify 
an expansion in powers of the shear rate for both these 
quantities about their equilibrium values 
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where due to the invariance of the shear direction only 
even powers in the shear rate 7 need to be considered. 

If we combine these expansions with the expression 
from classical theory one can easily derive expressions for 
AG*, the height of the nucleation barrier and N* , the size 
of the critical nucleus, both of which depend quadrati- 
cally on the shear rate. In Fig. |21 we show the results 
from our simulations where we extracted the height of the 
nucleation barrier for various pressures and shear rates. 
The dependence on the shear rate is confirmed by the 
parabolic fits. However, we caution the reader that this 
observation should not be considered as evidence that the 
shear rate can really be considered as a thermodynamic 
variable. In fact, in a recent study of the effect of shear on 
the location of the solid-liquid coexistence in a Lennard- 
Jones system, Butler and Harrowell found that no purely 
thermodynamic description of the effect of shear was pos- 
sible 27]. Shear directly affects the transport of particles 
from the solid to the liquid phase, and this effect is not 
thermodynamic. The expansion in Eq. Q is simply a 
way to represent the effect of shear as if it were purely 
thermodynamic. With this caveat in mind, we continue 
the remainder of the discussion in the language of clas- 
sical nucleation theory. We find that N* , the number of 
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FIG. 2: The height of the nucleation barrier /3AG* as function 
of the dimensionless shear rate 7/K^Z) for different pressures 
P. The solid lines are parabolic fits through the data. 

particles inside the critical cluster, also depends quadrat- 
ically on the applied shear rate. Using the classical nucle- 
ation theory expressions N* = (327r7^)/(3p||A/Ltp) and 
AG* = iV* I Ayu|/2, we can obtain the values of the second 
order coefficients in Eq. Q from a fit of the simulation 
data. The results are summarized in Table Q] We find a 



negative Cq implying a destabilization of the solid upon 
shear and a relatively small correction of the interfacial 
free energy. Both effects do not strongly depend on the 
pressure. Note, however, that the fits for A/x and 7 do 
not yield a good prediction for the shape of the nucle- 
ation barrier. The shape shows deviations from the one 
expected by classical nucleation theory, which is due to 
finite size effects of the cluster. 
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TABLE I: Numerical data for different pressures pP/n^ on the 
equilibrium barrier height AG'°''\ critical nucleus size A^''^''', 
and second order corrections to the free energy difference and 
interfacial free energy as obtained from the fitted simulation 
data. 

A bond order analysis shows that the structure of 
the nucleus is predominantly body-centered-cubic. Since 
small nuclei are in general neither spherical nor compact 
we have chosen to characterize their shape by the three 
principal moments of inertia. For a truly spherical nu- 
cleus these values would be identical, but since the shape 
of nucleus is fluctuating these moments are different. For 
relatively small clusters of 100 particles the ratio of the 
principle moments is roughly 6:10:12. As the nuclei grow 
larger, the differences between these moments of inertia 
get somewhat less. Surprisingly, the imposition of shear 
does not influence these ratios. This leads to the conclu- 
sion that although the size of the critical nucleus increases 
with shear, the overall shape is hardly influenced. This is 
different from the radial distribution functions we mea- 
sured in the liquid under shear. They become increas- 
ingly asymmetric for higher shear rates |23| . 

Knowledge of the eigenvectors of the inertia tensor al- 
lows us to determine its orientation. We find that the 
average orientation of the nucleus is weakly coupled to 
the direction of the applied shear. In particular, we find 
that the axis with the largest principal moment of iner- 
tia is, preferably in the gradient direction, in qualitative 
difference to a typical nearest neighbor particle cluster in 
a sheared fluid that prefers to be in the shear direction. 
The axis of the smallest principal moment of the nucleus 
tends to align with the vorticity direction. This align- 
ment becomes more pronounced with increasing nucleus 
size and with increasing shear rate. 

In Fig. 13 we show the orientation of the nucleus with 
respect to the shear direction. The tilt angle increases 
linearly with the applied shear rate 7 and only depends 
weakly on the osmotic pressure. In order to improve the 
statistical accuracy we have averaged over all cluster sizes 
between N = 100 and the critical nucleus size, N*. The 
inset of Fig.Olshows a schematic drawing of the preferred 
orientation of a nucleus. Note that the largest dimension 
of the nucleus (smallest principle moment) is preferably 



along the vorticity direction, i.e. perpendicular to the 
plane of drawing. Interestingly, a similar tilt occurs when 
vesicles with a flexible shape are exposed to a linear shear 
flow 29]. 
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FIG. 3: The tilt angle 9 of the principal moment of inertia 
with respect to the j/-axis. The inset shows a schematic rep- 
resentation of the preferred orientation of the nucleus with 
respect to the shear direction indicated by the arrows. 

In conclusion, we applied the combination of umbrella 
sampling and Brownian Dynamics simulation to the non- 
equilibrium problem of nucleation under shear, and found 
that shear suppresses nucleation and leads to a larger 
critical nucleus. These results can be described (but not 
yet understood) using a naive extension of classical nu- 
cleation theory. Most importantly, the present numeri- 
cal predictions can be tested experimentally, by studying 
the rate of homogeneous crystal nucleation in a homoge- 
neously sheared colloidal suspension. If nucleation were 
to be studied in Poiseuille flow as realized in a capil- 
lary viscometer [SOJ , rather than in homogeneous Couette 
flow, we should expect crystal nuclei to appear preferen- 
tially in the middle of the flow channel. 

We stress that the present flndings apply to the case 
where the fluid is only weakly sheared, i.e. when shear- 
induced ordering in the liquid phase is, presumably, 
unimportant. We also note that the present results in- 
dicate that, during sedimentation of crystal nuclei in an 
otherwise stagnant solution, local shear should decrease 
the rate of growth of the crystallites. There may even be 
conditions where the competition between mass gain due 
to crystal growth and mass loss due to shearing, leads 
to the selection of one particular crystallite radius. This 
phenomenon should also be experimentally observable. 

In this work, we ignored hydrodynamic interactions be- 
cause otherwise the computational cost would have been 
prohibitive. This assumption, while reasonable for dilute 
suspensions of charged colloids, is certainly not correct 
in general. Finally, our method can readily be applied 
to other dynamical simulation methods for rare events 
and meta-stable systems, such as crystal nucleation in 



oscillatory shear [3l| and heterogeneous nucleation near 
a system wall in a sheared suspension. 
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